#clean up
rm(list=ls())

#Set Working Directory: REVISE THIS TO WHERE YOU SAVE THE DATA.
setwd('/Volumes/MONOGAN/mediaCount/data/')

#Load Libraries
library(MASS)
library(TSA)
library(foreign)
library(tseries)
library(lattice)

#Load Patrick Brandt's code for a Poisson autoregressive model:
source('http://www.utdallas.edu/~pbrandt/code/pests.r')

########URA 2009, TABLE 3##############
ura<-read.csv('UraData.csv')

pdf("uraHist.pdf",family='Times',width=6,height=6,pointsize=10)
#postscript("uraHist.eps",family='ComputerModern',width=6,height=6,pointsize=10)
#histogram(~ura$usat,xlab="Homosexuality (USA Today)")
plot(table(ura$usat),type='h',axes=F,xlab="Homosexuality (USA Today)",ylab="Frequency",xlim=c(0,120))
axis(1,at=seq(0,120,by=10))
axis(2)
abline(h=0,col='gray60')
dev.off()

###Negative Binomial Regression###
ura.table.3.nb<-glm.nb(ura$usat~ura$romer+ura$lawr_long); summary(ura.table.3.nb)

#correlation between true and fitted
cor(ura$usat[-c(1:3)],ura.table.3.nb$fitted.values[-c(1:3)])
cor(ura$usat[-c(1:3)],ura.table.3.nb$fitted.values[-c(1:3)])^2

#RMSE
res<-ura$usat-ura.table.3.nb$fitted.values
sqrt(sum(res[-c(1:3)]^2))

#LB
Box.test(ura.table.3.nb$residuals,16,"Ljung-Box")

#Jarque-Bera
jarque.bera.test(ura.table.3.nb$residuals)

#long-run effects
100*(exp(ura.table.3.nb$coefficients)-1)

###Identification and PAR Model###
#standardized residuals from NB reg
ura.std<-(ura$usat-ura.table.3.nb$fitted.values)/sqrt(ura.table.3.nb$fitted.values)

##Self-Identify AR process##
acf(ura.std)
pacf(ura.std)
#Looks like an AR(3).

##PAR Model##
ura.table.3.par3<-Parp(ura$usat~ura$romer+ura$lawr_long, p=3); summary(ura.table.3.par3)
ura.par3.std<-(ura.table.3.par3$residuals)/sqrt(ura$usat-ura.table.3.par3$residuals)
acf(ura.par3.std,26)
pacf(ura.par3.std,26)
Box.test(ura.par3.std,26,'Ljung-Box')
Box.test(ura.par3.std,16,'Ljung-Box')

#correlation between true and fitted
par.fit<-ura$usat-ura.table.3.par3$residuals
cor(ura$usat[-c(1:3)],par.fit[-c(1:3)])
cor(ura$usat[-c(1:3)],par.fit[-c(1:3)])^2

#RMSE
sqrt(sum(ura.table.3.par3$residuals[-c(1:3)]^2))

#LB
Box.test(ura.par3.std,16,'Ljung-Box')
Box.test(ura.table.3.par3$residuals,16,"Ljung-Box")

#Jarque-Bera
jarque.bera.test(ura.table.3.par3$residuals)

parp.multipliers(ura.table.3.par3)

#long-run effects
parp.multipliers(ura.table.3.par3,xvec=matrix(c(1,0,1),1,ura.table.3.par3$k))
100*((19.0010520652343686265/par.fit[161]))

###Transfer Function###
ura.table.3.orig<-arimax(ura$usat, order=c(1,0,1), xtransf=cbind(ura$romer,ura$lawr_long), transfer=list(c(1,0),c(1,0))); ura.table.3.orig
ura.table.3.ar3<-arimax(ura$usat, order=c(3,0,0), xtransf=cbind(ura$romer,ura$lawr_long), transfer=list(c(1,0),c(1,0))); ura.table.3.ar3

#correlation between true and fitted
tf.fit<-ura$usat-ura.table.3.ar3$residuals
cor(ura$usat[-c(1:3)],tf.fit[-c(1:3)])
cor(ura$usat[-c(1:3)],tf.fit[-c(1:3)])^2

#RMSE
sqrt(sum(ura.table.3.ar3$residuals[-c(1:3)]^2))

#LB
Box.test(ura.table.3.ar3$residuals,16,"Ljung-Box")

#Jarque-Bera
jarque.bera.test(ura.table.3.ar3$residuals)

#long-run effects
24.250102616617930096/(1.23100067079911235091)
100*((19.699503982297436266/tf.fit[161]))

###OLS with 3 lags. (Koyck-like.)###
t.usat<-ts(ura$usat)
t.lawr_long<-ts(ura$lawr_long)
t.romer<-ts(ura$romer)
l.usat<-lag(t.usat,-1)
l2.usat<-lag(t.usat,-2)
l3.usat<-lag(t.usat,-3)
ura.2<-ts.union(t.usat,t.lawr_long,t.romer,l.usat,l2.usat,l3.usat)
ura.table.3.koyck<-lm(t.usat~l.usat+l2.usat+l3.usat+t.romer+t.lawr_long, data=ura.2); summary(ura.table.3.koyck)

#correlation between true and fitted
cor(ura$usat[-c(1:3)],ura.table.3.koyck$fitted.values)
cor(ura$usat[-c(1:3)],ura.table.3.koyck$fitted.values)^2

#RMSE
sqrt(sum(ura.table.3.koyck$residuals^2))

#LB
Box.test(ura.table.3.koyck$residuals,16,"Ljung-Box")

#Jarque-Bera
jarque.bera.test(ura.table.3.koyck$residuals)

#long-run effects
1.802063908602899556/(1-0.343820446592114193-0.241351662644831821-0.216150253522300051)
100*((9.0702906156428237949/ura.table.3.koyck$fitted.values[161]))

###Log-Normal Model###
log.usat.2<-ts(log(ura$usat+.1))
l.log.usat<-lag(log.usat.2,-1)
l2.log.usat<-lag(log.usat.2,-2)
l3.log.usat<-lag(log.usat.2,-3)
ura.3<-ts.union(log.usat.2,l.log.usat,l2.log.usat,l3.log.usat,t.lawr_long,t.romer)
ura.table.3.lognorm<-lm(log.usat.2~l.log.usat+l2.log.usat+l3.log.usat+t.romer+t.lawr_long, data=ura.3); summary(ura.table.3.lognorm)

#correlation between true and fitted
cor(ura$usat[-c(1:3)],exp(ura.table.3.lognorm$fitted.values))
cor(ura$usat[-c(1:3)],exp(ura.table.3.lognorm$fitted.values))^2

#RMSE
res<-ura$usat[-c(1:3)]-exp(ura.table.3.lognorm$fitted.values)
sqrt(sum(res^2))

#LB
Box.test(ura.table.3.lognorm$residuals,16,"Ljung-Box")

#Jarque-Bera
jarque.bera.test(ura.table.3.lognorm$residuals)

#long-run effects
100*(exp(0.048821133842861249/(1-0.386530225809196049-0.244032694930459815-0.162412511671727933))-1)

########INPUT: Romer v. Evans########
###Clarify-like simulations###
ura.3.sim.par<-mvrnorm(n=1000,mu=ura.table.3.par3$coefs[,1],Sigma=ura.table.3.par3$covar,empirical=TRUE)

ura.3.results.nb<-summary(ura.table.3.nb)
ura.3.sim.nb<-mvrnorm(n=1000,mu=ura.3.results.nb$coefficients[,1],Sigma=ura.3.results.nb$cov.unscaled,empirical=TRUE)

ura.3.results.koyck<-summary(ura.table.3.koyck)
ura.3.covar.koyck<-ura.3.results.koyck$cov.unscaled*(ura.3.results.koyck$sigma^2)
ura.3.sim.koyck<-mvrnorm(n=1000,mu=ura.3.results.koyck$coefficients[,1],Sigma=ura.3.covar.koyck,empirical=TRUE)

ura.3.results.lognorm<-summary(ura.table.3.lognorm)
ura.3.covar.lognorm<-ura.3.results.lognorm$cov.unscaled*(ura.3.results.lognorm$sigma^2)
ura.3.sim.lognorm<-mvrnorm(n=1000,mu=ura.3.results.lognorm$coefficients[,1],Sigma=ura.3.covar.lognorm,empirical=TRUE)

ura.3.sim.tf<-mvrnorm(n=1000,mu=ura.table.3.ar3$coef,Sigma=ura.table.3.ar3$var.coef,empirical=TRUE)

#apply(ura.3.sim.par,2,mean)
#apply(ura.3.sim.par,2,sd)
# 'usat' in April 1996=42. This is what all pre-intervention values are set to.

###Transfer function###
#Start with predictions from reported coefficients.
coef<-ura.table.3.ar3$coef[-c(1:3,5,7)]
romer.initial<-ura.table.3.ar3$coef[6]
romer.decay<-ura.table.3.ar3$coef[5]
lawr.initial<-ura.table.3.ar3$coef[8]
lawr.decay<-ura.table.3.ar3$coef[7]

#Romer onsets at wave 1. Waves -1-10
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0)
x1<-c(1,1,0)
y00<-coef%*%x00
y0<-coef%*%x0
y1<-coef%*%x1
y2<-coef%*%x2+romer.decay*romer.initial
y3<-coef%*%x3+(romer.decay^2)*romer.initial
y4<-coef%*%x4+(romer.decay^3)*romer.initial
y5<-coef%*%x5+(romer.decay^4)*romer.initial
y6<-coef%*%x6+(romer.decay^5)*romer.initial
y7<-coef%*%x7+(romer.decay^6)*romer.initial
y8<-coef%*%x8+(romer.decay^7)*romer.initial
y9<-coef%*%x9+(romer.decay^8)*romer.initial
y10<-coef%*%x10+(romer.decay^9)*romer.initial

#Plot Forecast
Time<-c(-1:10)
mle.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=mle.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
coef<-ura.3.sim.tf[i,-c(1:3,5,7)]
romer.initial<-ura.3.sim.tf[i,6]
romer.decay<-ura.3.sim.tf[i,5]
lawr.initial<-ura.3.sim.tf[i,8]
lawr.decay<-ura.3.sim.tf[i,7]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0)
x1<-c(1,1,0)
forecast[i,1]<-y00<-coef%*%x00
forecast[i,2]<-y0<-coef%*%x0
forecast[i,3]<-y1<-coef%*%x1
forecast[i,4]<-y2<-coef%*%x2+romer.decay*romer.initial
forecast[i,5]<-y3<-coef%*%x3+(romer.decay^2)*romer.initial
forecast[i,6]<-y4<-coef%*%x4+(romer.decay^3)*romer.initial
forecast[i,7]<-y5<-coef%*%x5+(romer.decay^4)*romer.initial
forecast[i,8]<-y6<-coef%*%x6+(romer.decay^5)*romer.initial
forecast[i,9]<-y7<-coef%*%x7+(romer.decay^6)*romer.initial
forecast[i,10]<-y8<-coef%*%x8+(romer.decay^7)*romer.initial
forecast[i,11]<-y9<-coef%*%x9+(romer.decay^8)*romer.initial
forecast[i,12]<-y10<-coef%*%x10+(romer.decay^9)*romer.initial
}

#Lower and Upper Bounds
mle.p.5<-apply(forecast, 2, quantile,.05)
mle.p.95<-apply(forecast, 2, quantile,.95)


###PAR model###
#Start with predictions from reported coefficients.
rho<-ura.table.3.par3$coefs[1:3,1]
weight<-1-sum(rho)
beta<-ura.table.3.par3$coefs[4:6,1]

#Romer onsets at wave 1. Waves -1-10
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0)
x1<-c(1,1,0)
lag00<-c(42,42,42)
y00<-weight*exp(beta%*%x00)+rho%*%lag00
lag0<-c(42,42,42)
y0<-weight*exp(beta%*%x0)+rho%*%lag0
lag1<-c(42,42,42)
y1<-weight*exp(beta%*%x1)+rho%*%lag1
lag2<-c(y1,42,42)
y2<-weight*exp(beta%*%x2)+rho%*%lag2
lag3<-c(y2,y1,42)
y3<-weight*exp(beta%*%x3)+rho%*%lag3
lag4<-c(y3,y2,y1)
y4<-weight*exp(beta%*%x4)+rho%*%lag4
lag5<-c(y4,y3,y2)
y5<-weight*exp(beta%*%x5)+rho%*%lag5
lag6<-c(y5,y4,y3)
y6<-weight*exp(beta%*%x6)+rho%*%lag6
lag7<-c(y6,y5,y4)
y7<-weight*exp(beta%*%x7)+rho%*%lag7
lag8<-c(y7,y6,y5)
y8<-weight*exp(beta%*%x8)+rho%*%lag8
lag9<-c(y8,y7,y6)
y9<-weight*exp(beta%*%x9)+rho%*%lag9
lag10<-c(y9,y8,y7)
y10<-weight*exp(beta%*%x10)+rho%*%lag10

#Plot Forecast
Time<-c(-1:10)
par.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=par.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
beta<-ura.3.sim.par[i,4:6]
rho<-ura.3.sim.par[i,1:3]
weight<-1-sum(rho)
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0)
x1<-c(1,1,0)
lag00<-c(42,42,42)
forecast[i,1]<-y00<-weight*exp(beta%*%x00)+rho%*%lag00
lag0<-c(42,42,42)
forecast[i,2]<-y0<-weight*exp(beta%*%x0)+rho%*%lag0
lag1<-c(42,42,42)
forecast[i,3]<-y1<-weight*exp(beta%*%x1)+rho%*%lag1
lag2<-c(y1,42,42)
forecast[i,4]<-y2<-weight*exp(beta%*%x2)+rho%*%lag2
lag3<-c(y2,y1,42)
forecast[i,5]<-y3<-weight*exp(beta%*%x3)+rho%*%lag3
lag4<-c(y3,y2,y1)
forecast[i,6]<-y4<-weight*exp(beta%*%x4)+rho%*%lag4
lag5<-c(y4,y3,y2)
forecast[i,7]<-y5<-weight*exp(beta%*%x5)+rho%*%lag5
lag6<-c(y5,y4,y3)
forecast[i,8]<-y6<-weight*exp(beta%*%x6)+rho%*%lag6
lag7<-c(y6,y5,y4)
forecast[i,9]<-y7<-weight*exp(beta%*%x7)+rho%*%lag7
lag8<-c(y7,y6,y5)
forecast[i,10]<-y8<-weight*exp(beta%*%x8)+rho%*%lag8
lag9<-c(y8,y7,y6)
forecast[i,11]<-y9<-weight*exp(beta%*%x9)+rho%*%lag9
lag10<-c(y9,y8,y7)
forecast[i,12]<-y10<-weight*exp(beta%*%x10)+rho%*%lag10
}

#Lower and Upper Bounds
par.p.5<-apply(forecast, 2, quantile,.05)
par.p.95<-apply(forecast, 2, quantile,.95)


###Negative binomial regression###
#Start with predictions from reported coefficients.
coef<-ura.3.results.nb$coefficients[,1]

#Romer onsets at wave 1. Waves -1-10
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0)
x1<-c(1,1,0)
y00<-exp(coef%*%x00)
y0<-exp(coef%*%x0)
y1<-exp(coef%*%x1)
y2<-exp(coef%*%x2)
y3<-exp(coef%*%x3)
y4<-exp(coef%*%x4)
y5<-exp(coef%*%x5)
y6<-exp(coef%*%x6)
y7<-exp(coef%*%x7)
y8<-exp(coef%*%x8)
y9<-exp(coef%*%x9)
y10<-exp(coef%*%x10)

#Plot Forecast
Time<-c(-1:10)
nb.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=nb.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
coef<-ura.3.sim.nb[i,]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0)
x1<-c(1,1,0)
forecast[i,1]<-y00<-exp(coef%*%x00)
forecast[i,2]<-y0<-exp(coef%*%x0)
forecast[i,3]<-y1<-exp(coef%*%x1)
forecast[i,4]<-y2<-exp(coef%*%x2)
forecast[i,5]<-y3<-exp(coef%*%x3)
forecast[i,6]<-y4<-exp(coef%*%x4)
forecast[i,7]<-y5<-exp(coef%*%x5)
forecast[i,8]<-y6<-exp(coef%*%x6)
forecast[i,9]<-y7<-exp(coef%*%x7)
forecast[i,10]<-y8<-exp(coef%*%x8)
forecast[i,11]<-y9<-exp(coef%*%x9)
forecast[i,12]<-y10<-exp(coef%*%x10)
}

#Lower and Upper Bounds
nb.p.5<-apply(forecast, 2, quantile,.05)
nb.p.95<-apply(forecast, 2, quantile,.95)

###OLS with 3 lags. (Koyck-like.)###
#Start with predictions from reported coefficients.
coef<-ura.3.results.koyck$coefficients[,1]

#Romer onsets at wave 1. Waves -1-10
x00<-c(1,42,42,42,0,0)
y00<-coef%*%x00
x0<-c(1,42,42,42,0,0)
y0<-coef%*%x0
x1<-c(1,42,42,42,1,0)
y1<-coef%*%x1
x2<-c(1,y1,42,42,0,0)
y2<-coef%*%x2
x3<-c(1,y2,y1,42,0,0)
y3<-coef%*%x3
x4<-c(1,y3,y2,y1,0,0)
y4<-coef%*%x4
x5<-c(1,y4,y3,y2,0,0)
y5<-coef%*%x5
x6<-c(1,y5,y4,y3,0,0)
y6<-coef%*%x6
x7<-c(1,y6,y5,y4,0,0)
y7<-coef%*%x7
x8<-c(1,y7,y6,y5,0,0)
y8<-coef%*%x8
x9<-c(1,y8,y7,y6,0,0)
y9<-coef%*%x9
x10<-c(1,y9,y8,y8,0,0)
y10<-coef%*%x10

#Plot Forecast
Time<-c(-1:10)
ols.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=ols.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
coef<-ura.3.sim.koyck[i,]
x00<-c(1,42,42,42,0,0)
forecast[i,1]<-y00<-coef%*%x00
x0<-c(1,42,42,42,0,0)
forecast[i,2]<-y0<-coef%*%x0
x1<-c(1,42,42,42,1,0)
forecast[i,3]<-y1<-coef%*%x1
x2<-c(1,y1,42,42,0,0)
forecast[i,4]<-y2<-coef%*%x2
x3<-c(1,y2,y1,42,0,0)
forecast[i,5]<-y3<-coef%*%x3
x4<-c(1,y3,y2,y1,0,0)
forecast[i,6]<-y4<-coef%*%x4
x5<-c(1,y4,y3,y2,0,0)
forecast[i,7]<-y5<-coef%*%x5
x6<-c(1,y5,y4,y3,0,0)
forecast[i,8]<-y6<-coef%*%x6
x7<-c(1,y6,y5,y4,0,0)
forecast[i,9]<-y7<-coef%*%x7
x8<-c(1,y7,y6,y5,0,0)
forecast[i,10]<-y8<-coef%*%x8
x9<-c(1,y8,y7,y6,0,0)
forecast[i,11]<-y9<-coef%*%x9
x10<-c(1,y9,y8,y8,0,0)
forecast[i,12]<-y10<-coef%*%x10
}

#Lower and Upper Bounds
ols.p.5<-apply(forecast, 2, quantile,.05)
ols.p.95<-apply(forecast, 2, quantile,.95)

###Log-Normal Model###
#Start with predictions from reported coefficients.
coef<-ura.3.results.lognorm$coefficients[-c(2:4),1]
rho1<-ura.3.results.lognorm$coefficients[2,1]
rho2<-ura.3.results.lognorm$coefficients[3,1]
rho3<-ura.3.results.lognorm$coefficients[4,1]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0)
x1<-c(1,1,0)
y00<-exp(coef%*%x00)*(42^rho1)*(42^rho2)*(42^rho3)
y0<-exp(coef%*%x0)*(42^rho1)*(42^rho2)*(42^rho3)
y1<-exp(coef%*%x1)*(42^rho1)*(42^rho2)*(42^rho3)
y2<-exp(coef%*%x2)*(y1^rho1)*(42^rho2)*(42^rho3)
y3<-exp(coef%*%x3)*(y2^rho1)*(y1^rho2)*(42^rho3)
y4<-exp(coef%*%x4)*(y3^rho1)*(y2^rho2)*(y1^rho3)
y5<-exp(coef%*%x5)*(y4^rho1)*(y3^rho2)*(y2^rho3)
y6<-exp(coef%*%x6)*(y5^rho1)*(y4^rho2)*(y3^rho3)
y7<-exp(coef%*%x7)*(y6^rho1)*(y5^rho2)*(y4^rho3)
y8<-exp(coef%*%x8)*(y7^rho1)*(y6^rho2)*(y5^rho3)
y9<-exp(coef%*%x9)*(y8^rho1)*(y7^rho2)*(y6^rho3)
y10<-exp(coef%*%x10)*(y9^rho1)*(y8^rho2)*(y7^rho3)

#Plot Forecast
Time<-c(-1:10)
lognorm.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=lognorm.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
coef<-ura.3.sim.lognorm[i,-c(2:4)]
rho1<-ura.3.sim.lognorm[i,2]
rho2<-ura.3.sim.lognorm[i,3]
rho3<-ura.3.sim.lognorm[i,4]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0)
x1<-c(1,1,0)
forecast[i,1]<-y00<-exp(coef%*%x00)*(42^rho1)*(42^rho2)*(42^rho3)
forecast[i,2]<-y0<-exp(coef%*%x0)*(42^rho1)*(42^rho2)*(42^rho3)
forecast[i,3]<-y1<-exp(coef%*%x1)*(42^rho1)*(42^rho2)*(42^rho3)
forecast[i,4]<-y2<-exp(coef%*%x2)*(y1^rho1)*(42^rho2)*(42^rho3)
forecast[i,5]<-y3<-exp(coef%*%x3)*(y2^rho1)*(y1^rho2)*(42^rho3)
forecast[i,6]<-y4<-exp(coef%*%x4)*(y3^rho1)*(y2^rho2)*(y1^rho3)
forecast[i,7]<-y5<-exp(coef%*%x5)*(y4^rho1)*(y3^rho2)*(y2^rho3)
forecast[i,8]<-y6<-exp(coef%*%x6)*(y5^rho1)*(y4^rho2)*(y3^rho3)
forecast[i,9]<-y7<-exp(coef%*%x7)*(y6^rho1)*(y5^rho2)*(y4^rho3)
forecast[i,10]<-y8<-exp(coef%*%x8)*(y7^rho1)*(y6^rho2)*(y5^rho3)
forecast[i,11]<-y9<-exp(coef%*%x9)*(y8^rho1)*(y7^rho2)*(y6^rho3)
forecast[i,12]<-y10<-exp(coef%*%x10)*(y9^rho1)*(y8^rho2)*(y7^rho3)
}

#Lower and Upper Bounds
lognorm.p.5<-apply(forecast, 2, quantile,.05)
lognorm.p.95<-apply(forecast, 2, quantile,.95)

###Plot All Forecasts With Error Bands###
pdf("uraRomer.pdf",family='Times',width=6,height=6,pointsize=10)
#postscript("uraRomer.eps",family='ComputerModern',width=6,height=6,pointsize=10)
Time<-c(-1:10)

#points(y=ols.Forecast,x=Time+.1,pch=23,col='blue')
plot(y=ols.Forecast,x=Time-.26,pch=23,ylim=c(30,140),xlim=c(-1,11),xlab="Time",ylab="Predicted Count",col='blue', axes=F)
arrows(x0=Time-.26,y0=ols.p.5,x1=Time-.26,y1=ols.p.95, code=3, angle=90, length=.0,lty=2,col='blue')
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('April 96','June 96','Aug 96','Oct 96','Dec 96','Feb 97'))
box()

#plot(y=nb.Forecast,x=Time-.1,pch=20,ylim=c(0,250),xlim=c(-1,11),xlab="Time",ylab="Predicted Count",col='red', axes=F)
points(y=nb.Forecast,x=Time-.13,pch=20,col='red')
arrows(x0=Time-.13,y0=nb.p.5,x1=Time-.13,y1=nb.p.95, code=3, angle=90, length=.0,lty=3,col='red')

points(y=lognorm.Forecast,x=Time,pch=15,col='purple')
arrows(x0=Time,y0=lognorm.p.5,x1=Time,y1=lognorm.p.95, code=3, angle=90, length=.0,lty=5,col='purple')

points(y=mle.Forecast,x=Time+.13,pch=22,col='forestgreen')
arrows(x0=Time+.13,y0=mle.p.5,x1=Time+.13,y1=mle.p.95, code=3, angle=90, length=.0,lty=4,col='forestgreen')

points(y=par.Forecast,x=Time+.26,pch=24)
arrows(x0=Time+.26,y0=par.p.5,x1=Time+.26,y1=par.p.95, code=3, angle=90, length=.0,lty=1)

legend(x=2,y=135,legend=c('Koyck','Negative Binomial','Log-normal','Transfer Function','PAR'), lty=c(2,3,5,4,1), pch=c(23,20,15,22,24), col=c('blue','red','purple','forestgreen','black'))
dev.off()

